Analytic Model of the Universal Structure of Turbulent Boundary Layers 
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Turbulent boundary layers exhibit a universal structure which nevertheless is rather complex, 
being composed of a viscous sub-layer, a buffer zone, and a turbulent log-law region. In this letter we 
present a simple analytic model of turbulent boundary layers which culminates in explicit formulae 
for the profiles of the mean velocity, the kinetic energy and the Reynolds stress as a function of the 
distance from the wall. The resulting profiles are in close quantitative agreement with measurements 
over the entire structure of the boundary layer, without any need of re-fitting in the different zones. 
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Introudction: Theoretical physicists tend to con- 
sider turbulence in the context of the idealized model of 
isotropic and homogenous fluid flows at large Reynolds 
numbers. In part this is due to the apparent existence of 
universal, anomalous scaling exponents which character- 
ize correlation and structure functions in fully developed 
turbulent flows. It is also easier to disregard the effects 
of walls which introduce essential anisotropics and inho- 
mogeneities. Needless to say, all realistic turbulent flows 
are neither homogeneous nor isotropic. A problem of 
extreme interest for both technological applications and 
from the point of view of basic science is "wall-bounded" 
turbulence, with the theoretical model of a flat infinite 
wall playing a key role. This problem presents also fas- 
cinating universal features, see Fig.l, but traditionally it 
was more popular in the engineering rather than in the 
physics community. There are fascinating open problems 
in wall-bounded turbulence. The present Letter attempts 
at finding a simple model that affords an analytic calcu- 
lation of the universal profiles of the mean velocity, the 
turbulent energy and the Reynolds stress as a function 
of the distance from the wall. 

The theory that we construct begins with the equa- 
tions of fluid mechanics and focuses on the momentum 
and energy fluxes using the conservation laws for these 
quantities as a guidance for developing an appropriate 
model. In gross substance this approach is not new, and 
indeed a number of ingredients are borrowed from the 
literature. The model that we end up with is however 
improved compared to previous results in the sense that 
it provides analytic predictions for the above-mentioned 
profiles in the entire boundary layer, without re-fitting 
in the different zones. It should be stressed at this point 
that one cannot expect a universal model to apply for 
all turbulent boundary values problem. For instance, 
the understanding of drag reduction by additives calls 
for a slightly modified model that stresses the effects of 
the additives. The present contribution offers what we 
consider a simplest model that is constructed to best de- 
scribe Newtonian wall-bounded turbulence with enough 
richness to capture all the essential universal profiles of 
the quantities of interest. One advantage of the model 
that will be demonstrated in future publications is that it 
can be naturally generalized to describe stratified turbu- 
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FIG. 1: A sketch of the characteristic regions appearing in a 
turbulent boundary layer. The continuous line represents the 
mean velocity profile V'^{y'^) [in the wall units, Eq. (4)] taken 
from direct numerical simulation at Re^ = 2320, [1]. The 
"viscous region" ends when the linear law — begins 
to deviate from the continuous line at y^ ~ 6. The "log-law 
region" ends approximately at 1/6 of the channel hight, i.e. 
y'^— Rbt/O. At this point and deeper towards the channel 
cenerline the dimensionless momentum flux to the wall, = 
1 — J/"*"/ Rer (see Eq. (24)), deviates appreciable from unity. 



lent boundary layers caused by, say, temperature gradi- 
ents, heavy particles, etc, with applications to interesting 
problems like sand storms in deserts, snow fall in a windy 
days or water flow in a heavily silt-laden rivers. 

Equations and definitions: The starting point are 
the Navier-Stokes equations for an incompressible fluid 
velocity U{r,t), 

dU Vn 

^ + {U-V)U = -^ + iyAU, V-U = 0, (1) 
ot p 

where p is the fluid density, p — p{r, t) ~ the pressure 
and v is the kinematic viscosity. We follow the standard 
strategy of Reynolds considering the velocity as a sum of 
its average (over time) and a fluctuating part: 

U{r,t)^V{r)+u{r,t) , V{r) = {U{r,t)) . (2) 
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We also introduce the conventional viscous scale £r and 
friction velocity Ur 

Ur = VPiO)/P, tr = ylUr. (3) 

where P(0) is the wall shear stress the flux of the me- 
chanical momentum at the wall. This quantity depends 
on the type of turbulent boundary layer (pressure driven 

or driven by the upper boundary, etc). The velocity and 
the distance from the wall are then measured in wall units 

U+= U/ur , y+= y/£r , P+= p/{pul) , V+= 4V. (4) 

In wall units the Navier-Stokes equation is now dimen- 
sionless, and the stationary version takes the form 

(C/+ • V+) U+ = - V+P+ + A+U+ . (5) 

The averaged equation for the mean velocity looks like 

(y+ ■ V+) V+ = A+V+-{u+-V+u+)-V+{p+) . (6) 

Besides the mean velocity one needs to consider cor- 
relation functions. It turns out that important features 
of wall bounded turbulence, like the mean velocity pro- 
file, thermal conductivity, turbulent transport of matter, 
etc., are determined by velocity fluctuations on relatively 
large scales, and the statistics of the latter do not devi- 
ate too much from Gaussian statistics. Thus an economic 
description of these features can be reached on a level of 
second order correlation functions. Therefore, in addi- 
tion to the mean velocity profile, we will be interested in 
the detailed description of two additional quantities, the 
turbulent energy and the Reynolds stress tensor: 

K+ ^ (|n+|2)/2 , W± ^ {utu+) . (7) 

Balance equations for the Reynolds stress: Sub- 
tracting Eq. (6) from Eq. (5), multiplying by a fluctuat- 
ing velocity component and averaging, results in the 
equation satisfied by the Reynolds stress: 

{V+ . V+) W^^ = P++ 7^+ - e± + d,T+, , (8) 

where the tensors of the energy production, P^^, of the 
pressure-rate-of-strain, TZ^- and of the Reynolds-stress 
dissipation ej-, are 

Pt,^-Wid,V+ + W+,dkVt. (9) 
K ^ (P^ (.9,u+ + d,u+) ) , e+= 2{dkut dkUp , 

and p = p — (p) denotes the pressure fluctuations. The 
last term T.^^. presents spatial energy fluxes. We will 
neglect it throughout the turbulent boundary layer; it is 
indeed small in the log-layer, but comparable to other 
terms in the buffer and viscous sub-layers. The model 
will be constructed such as to compensate for this neglect 
in those regions where the term is significant. The bonus 



of neglecting this term is enormous since this keeps the 
theory local, without partial derivatives. 

The modeling of the various terms appearing in Eq. (8) 
has attracted considerable attention over the years, and 
we only briefly summarize how this is done. The Pois- 
son's equation for the fluctuating pressure follows from 
the equation of the fluctuating part of the velocity field, 
u, which is obtained by subtracting Eq. (6) from (5): 

A+P+ = -V+V+ {utu+ - (utul) + V+uj + V+ut) . 

(10a) 

The homogeneous solution of this equation is responsible 
for sound, a phenomenon of very little consequence for 
turbulent dynamics at low Mach numbers. The inhomo- 
geneous solution includes two parts, p~^ = p+j +Pu- 

pL <x utu+ - {u+u+} , p+ (X V^u+ + V+uf . (10b) 

Correspondingly the pressure-rate-of-strain tensor (9) 
consist of two terms: 

7^+ = i?-+ + i^-+. (iia) 

The first of these is known as the "Return-to-Isotropy" 
tensor, Rfj^, that depends on the tripple- velocity corre- 
lator (uiUjUk). Its evaluation in terms of the objects of 
the theory calls for a closure, and following time-honored 
tradition [2] we adopt for it the simple Rota form 

-7,, (3 W+ - W+) , W+ = Tr{W+} , (lib) 

in which 7^1 is some characteristic nonlinear frequency 
that will be specified later. The tensor Rif^ is traceless 
and therefore the frequency 7^1 must be the same for 
all diagonal components of R^'~^ ■ There is no reason 
however to assume that off-diagonal terms have the same 
frequency. Therefore, following [3], we assert that 

R^'+ ^ -3^n.Wij , i^j, (11c) 

with, generally speaking, 7^ 7^ 7^1. 

The traceless "Isotropization-of-Production" tensor, 
Rij^ has a structure that is very similar to the pro- 
duction tensor V^, Eq. (9), and thus traditionally it is 
modeled in terms of V^^ [2] : 

RIJ+ ^ -C:p (3P+ - P+Si,) , V+ ^ Tt{V+} . (lid) 

The dissipation tensor efj is estimated differently far 
from the wall and near it. Far from the wall and for large 
Reynolds numbers the turbulent flow can be considered 
approximately isotropic. Therefore, the tensor should 
be approximately diagonal, 

e+=j+W+Sij. (12a) 

Under stationary conditions the rate of turbulent kinetic 
energy dissipation is equal to the energy input at the 
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outer scale, estimated as (uiUjUk) j I where the outer 
scale of turbulence I is estimated as the distance to the 
wall y. Therefore, the natural estimate of 7+ involves 
the tripple-velocity correlator, 



(12b) 



Similarly, we can estimate the Return-to-Isotropy fre- 
quencies 7+ and 7+ in Eqs. (lib) and (11c). Having 
in mind that the precise structure (tensorial contraction, 
etc.) of the equation for 7+ is different from that of the 
equations for 7+ and -y^^, we should involve different nu- 
merical prefactors: 



i,=b^,Vw^/y+, ^+=b^,Vw^/y+. (13) 
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Close to the wall, in the viscous sub-layer t/+ < 30, the 
estimates change due to the direct viscous contribution 
of the largest eddies at that distance, whose size is of the 
order of the distance y+ itself. For these eddies we can 
estimate the operator as (a/y)"^, where a is a new 
fitting constant. For simplicity we account for this con- 
tribution only in off-diagonal terms, which do not include 
the nonlinear part (see Eq. (12a)). In this way we will 
have: 



£+. = r+ 



i 7^ J. 



(14a) 



with r+ = (a/y)'^. Here we have to recall that we ne- 
glected the spatial energy transfer term, which plays an 
important role in the viscous sub-layer. Since it has an 
opposite sign to the viscous dissipation, we can take its 
influence into account by suppressing the direct viscous 
dissipation by a a function of W/W^, (here is the 
asymptotical value of W in the lag-low region): 



r+ = {a/y+Y JW+/W+ 



(14b) 



Our choice of a square-root function is dictated by the 
simplicity of the analytical treatment of the resulting al- 
gebraic model. A-posteriori, the implicit accounting for 
the energy flux in this particular way is supported by the 
good agreement of the model prediction for the energy 
profile in both the viscous and the buffer layers with the 
DNS data without additional fitting parameters, see the 
insert in Fig. 3. 

Plane geometry and the balance of momentum: 
For plane geometry the mean velocity is oriented in the 
(stream-wise) x direction and depends only on the verti- 
cal (wall-normal) coordinate y: V = y(y)x. For such 
flows all the averages are functions of y"*" only. An 
interesting special example is a channel flow of height 
2L and infinite extent in the span-wise direction. Due 
to the symmetry in the span- wise direction z — > — z, 
= = 0. From Eq. (6) for V (integrated over 2/+) 



one gets the exact balance equation for the mechanical- 
momentum 

P+{y+)= S+-W+y, where (15) 

The flat geometry also simplifies the production term de- 
fined in Eq. (9): 



(16) 



Final set of equations: Substituting everything into 
Eq. (8) we get the following set of model equations: 

37iVt^.t = {it -l+)W+~2{l^2 C:p) S+W+y, (17a) 
^lt.W^y = (in. - 7+) W+ - 2 a,S+W+ , (17b) 
= (7+ - 7+) W+ - 2 a,S+W+ , (17c) 
= (r++ 37+) + (1 - 3 Cp) S+W+ . (17d) 

Summing up Eqs. (17a), (17b) and (17c) [or, equivalently, 
taking the trace of Eq. (8)] one gets 



3^+W+ = -2S+W+ 



(18) 



Choice of the log- law parameters b, 6ri and 6ri: 
At this point we can fit three 6-parameters responsible 
for the solution in the asymptotical log-law region 30 < 
2/+ < Re-r/S for sufficiently large Reynolds numbers. In 
this region r+ can be neglected, and solution of Eqs. (17) 
takes the form 



yy 
W+ = 

xy 



(l-2ap)6 + 6,: , 
6 + 36hi 


(19a) 


b + 3bni * 


(19b) 


-bK{W+f' = -l, 


(19c) 


— — , where 

Ky+ 


(19d) 


/2(1-3C:p)(C:p5 + 6„,) 

V 3W^bKi{b + 3b„,) 


(19e) 



To determine 6-parameters, we used the following data: 

1. The numerical values of = 7.8 and von-Karman 
constant k — 0.405 can be taken from the direct numer- 
ical simulations (DNS) [4]; 

2. The detailed analysis of experimental, DNS and 
large-eddy-simulation data, made in Ref. [3], yields the 
conclusion that with a good accuracy one can take 



W+=2W-j-=2W^, = W+/2; 



(20) 



3. The suggested in the literature (see, e.g. [2]) value 
of Cip is ^ . Our analysis showed that all profiles are very 
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insensitive to a particular choice of Cip around value 0.2. 
Therefore for simplicity we take 

C:p = 1/5 . (21a) 
All this knowledge enables us to find b, b^j and b^i: 

^ - -- (21b) 

(21c) 
(21d) 

(21e) 



w 0.075, 

4(1-3C:p)6w0.121, 
(l-3C,p) 



12 k 



0.230 



For future purpose we introduce a parameter 



1 - 3Cip 



12 k 



= 0.575 



General solution: With the chosen parameters (21) 
the formal solution of the system (15), (17) in the entire 
turbulent boundary layer is: 



xy 



2W^ =AW^ 

XX yy 



(22a) 



3 7+ 
2S+ 



^+'=15 7+(r+ + 3 7+). (22b) 



The last equation comes from the solvability condition 
for the system of Eqs. (17): Det = 0. Introduce: 



Then 



3 &m vt 



W. 



xy 




s+= 



3v^ 



V2bbr, 



and Eq. (15) transforms into: 



?;+ +12brv+/y+ - P+ J8br/b = 



(23) 



(24) 



(25) 



This is just a quadratic equation for w+ = with a 

unique positive solution: 
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6b 



(26) 



Comparison of the model and simulations: 

Clearly, a model with only 4 (or 5 if Cip is counted in) 
fit parameters cannot fit perfectly the profiles of all the 
physical quantities that can be measured. Therefore the 
actual value of the only remaining parameter a should 
be determined with a choice of the characteristics of tur- 
bulent boundary layers that we desire to describe best. 
Foremost in any modeling should be the mean velocity 
profile y which is of crucial importance in a wide variety 
of transport phenomena. Therefore we chose the value 




FIG. 2: Mean velocity profiles V {y ): The dotted line re- 
produces the results of direct numerical simulations [4] for 
Rer = 2000, the solid line is the analytical prediction of the 
model with a = 7.1 and the values (21) of h, bm and 6^. 



of a from the best fit of V^{y^) in the (quasi) straight 
logarithmic region 30 < ?/+ < 350: 



a = 7.1 



(27) 



The resulting mean velocity profile V^{y^) = 

Jo S+i£,)d^, in which S+{^) is given by (24), is shown 
in Fig. 2 by a solid line for Re^ = 2000. The dotted fine 
represents data taken from direct numerical simulations 
[4], for the same Re,-. There is no significant difference 
between these plots in the viscous sublayer, buffer and 
outer layers, where j/+ < 800 i.e. in about 40% of the 
channel half- width = Rct-. This robustness of the 
mean velocity profile is a consequence of the fact that 
y+(y+) is an integral of the mean shear 5+ which is de- 
scribed very well both in the viscous and the outer layers. 

Notice that our model does not describe the upward 
deviation from the log-low which is observed near the 
mid-channel (of about a few units in V^). We consider 
this minor disagreement as an acceptable price for the 
simplicity of the model which neglects the spatial energy 
transport term toward the centerline of the channel. This 
transport is the only reason for some turbulent activity 
near the centerline where both the Reynolds stress Wxy 
and S vanish due to symmetry. Just at the center line 
the source term in our energy equation, —2SWxy, is zero, 
and the missing energy transport term is felt. 

The plots in Fig. 2 have a reasonably straight loga- 
rithmic region from « 30 to y~^ « 350. On the other 
hand, the Reynolds stress profile at the same Re^ = 2000 
shown in Fig. 3, has no flat region at all. Such a 
flat region is expected in the true asymptotic regime of 
RCr oo, where W'^ = — 1. Therefore if one plots the 
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FIG. 3: The trace of the Reynolds-stress tensor (twice the 
total kinetic energy densiy): The dotted blue line reproduces 
the results of direct numerical simulations [4] at Rer = 2000, 
the solid red line is analytical prediction with for a = 7.1. 
The insert shows the buffer layer behavior in more detail. 
Notice, that there is no plateau in these plots, meaning that 
these values of Rer are not large enough to have a true scale- 
invariant log-law region. 

model profiles at different Re,- and fits them by log- 
linear profiles 

V+iy+) = K-Hn{y+)+B, (28) 

one can get a Rer-dependence of the "effective" inter- 
cept B =J> B(Iler) in the von-Karman log-law (28). We 
think that this explains why measured values of the in- 
tercept depend on the Reynolds number and on the flow 
geometry (channel vs. pipe): both in direct numerical 
simulations and in physical experiments one usually does 
not reach high enough values of Re,- . 

At this point a and all three b -parameters are already 
chosen from the information about the mean velocity pro- 
file and the ^"'"-independent values of the Reynolds stress 
tensor in the lag-law region in the limit Re,- — > oo. In 
Figs. 3 and 4 we show the predictions of the model for 
the kinetic energy profile W'^, and the Reynolds stress 
profile (red soHd lines). The result pertain to the entire 
turbulent boundary layer for ReT-=2000, without any ad- 
ditional fit parameters. In the same figures we show the 
results of direct numerical simulation [4] (blue dotted 
lines). The excellent agreement in the entire turbulent 
boundary layer is a strong indication that the proposed 
model summarizes the important physical ingredients of 
the problem. 

Summary: In summary, we presented an analytic 
model of the physics of wall-bounded turbulence in a 
Newtonian fluid, based entirely on the balance of energy 
and momentum fluxes with the production and the dissi- 
pation. The model has one blatant simplification which 
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FIG. 4: The Reynolds stress |W+,| : The dotted blue line 
represents results from direct numerical simulations [4] for 
Re^ — 2000, the solid red line represents the prediction of 
the analytical model with a = 7.1. There is no pronounced 
log- law region where a plateau |Wa^^| = 1 is expected. In- 
stead reaches only the value 0.85 around 50. This 
means that in this region the viscous transport is still essen- 
tial, but the total momentum flux is still below its (dimen- 
sionless) maximum value: max(7-'"'") = 1. 

is the neglect of the spatial energy fluxes. The gain is 
enormous - we get a local model that can be solved an- 
alytically to find the profiles of the mean velocity, the 
turbulent fluctuations and the Reyonlds stress as a func- 
tion of the distance from the wall, with all the boundary 
layer represented without re-fitting in the various regions. 
The reason of the success of this simple model is that we 
have learned how to compensate the neglect of the spatial 
flux in the buffer and the viscous region by a decrease in 
the dissipative terms. In future work we will demonstrate 
the utility of this simple model in a variety of important 
wall-bounded flows. 

Acknowledements: This research is supported by 
the US-Israel Binational Science Foundation. 



[1] K. Iwamoto, N. Kasagi, and Y. Suzuki, Direct numerical 

simulation of turbulent channel flow Re^ — 2320, Proc. 

6th Symp. Smart Control of Turbulence Tokyo, March 6- 

9, 2005, pp. 327-333. 
[2] S.B. Pope, Turbulent Flows, 1st ed. (Cambridge University 

Press, 2000). 

[3] V.S. L'vov, A. Pomyalov, I. Procaccia and S.S. Zilitinke- 
vich, Phys. Rev. E., 73, 016303 (2006). 

[4] S. Hoyas, J. Jimenez, Scaling of velocity fluctu- 
ations in turbulent channels up to Re^ = 2000, 
Phys. of Fluids, 18, 011702 (2006); DNS data at 
htt p : / /t orroj a. dmt.upm.es /ftp/ channels / . 



